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Abstract 


The present study is aimed at optimization of the geometric parameters of the proton exchange membrane (PEM) fuel cells through numerical 
simulation. The approach is developed by integrating a direct problem solver with an optimizer. A commercial computational fluid dynamics code 
is used as the direct problem solver, which is used to simulate the three-dimensional mass, momentum and species transport phenomena as well 
as the electron- and proton-transfer process taking place in a PEMFC. On the other hand, the simplified conjugate-gradient method (SCGM) is 
employed to build the optimizer, which is combined with the direct problem solver in order to seek the optimal geometric parameters, including, 
for example, the gas channel width fraction, the gas channel height and the thickness of the gas diffusion layer. It is found that the present approach 
can be applied to determine the optimal set of geometric parameters, and the search process is robust and always leads to a unique final solution 


regardless of the initial guess. 
© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Numerical methods have been employed to study the physi- 
cal phenomenon of interest in the fuel cells since these methods 
require relatively lower cost for analysis than the experimental 
ones. Furthermore, they in general provide a sufficient amount of 
data not easily obtainable with measurements. Thus, important 
natures of the flow and the thermal fields might then be observed 
and discussed extensively by numerical studies. To name a few, 
Wang and Savinell [1] performed a simulation of the effects of 
porosity and thickness of the catalyst layer, Pt loading and CO 
poisoning on fuel cells. Yi and Nguyen [2] used the numeri- 
cal methods to solve a two-dimensional single-phase PEMFC 
model with interdigitated flow channels so as to evaluate the 
effects of inlet and exit pressures, gas diffusion layer thickness 
and carbon plate width on the performance of PEMFC. Gurau 
et al. [3] developed a half-cell model to deal with the transport 
phenomena at the cathode of the PEMFC. The research explored 
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the effects of the porosity of the gas diffusion layer on the per- 
formance of fuel cells. A transient model has been developed by 
Umetal. [4] to simulate proton exchange membrane fuel cells by 
accounting simultaneously for electrochemical kinetics, current 
distribution, hydrodynamics and multicomponent transport. A 
two-dimensional cross-the-channel model was applied by Sun 
et al. [5] to investigate the influence of gas diffusion layer (GDL) 
property and flow-field geometry on the local reaction rate in the 
PEMFC cathode catalyst layer. 

As to the three-dimensional models, Hontanon et al. [6] 
adopted a commercial computational fluid dynamics code (FLU- 
ENT) to predict the effects of a group of geometrical parameters 
on the polarization performance of the fuel cells. Berning and 
Djilali [7] used a three-dimensional physical model to eval- 
uate the performance of the fuel cells at different operating 
conditions. Their research showed that an increase in the fuel 
cell temperature could raise the cell exchange current density, 
ionic conductivity of proton exchange membrane, gas diffu- 
sivity, and the operating voltage and then lead to an increase 
in the cell performance. A detailed review regarding the CFD 
modeling and related experimental validation has been given by 
Wang [8]. 
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Nomenclature 

a water activity 

D mass diffusivity (cm? s7!) 

e design parameter 

F Faraday constant, 96,487 C mol! 

h height of gas channel (m) 

i local electric current density (A m~?) 

i local ionic current density (A m7?) 

I electric current density (A m~?) 

j transfer current density (A m7?) 

J objective function 

k permeability (m7) 

Ic gas channel width (m) 

L gas channel length (m) 

lp width of a module (m) 

M molecular weight (kg kmol~!) 

n coordinate in the direction normal to the surface 

P pressure (Pa) 

P power density (W m~?) 

r concentration parameter 

R universal gas constant, 8.314 J mol~! K7! 

S source term in species equation 

tGDL thickness of gas diffusion layer (m) 

tCat thickness of catalyst layer (m) 

tm thickness of membrane (m) 

tp thickness of carbon plate (m) 

T temperature (K) 

U gas velocity vector (ms~!) 

V inlet gas velocity (m s7!) 

V cell voltage (V) 

X mole fraction 

x,y,z Cartesian coordinates 

Y mass fraction 

Greek symbols 

I'scat ionic conductivity of solid catalyst particle 
(Q7! m7!) 

In ionic conductivity of membrane (Q-! m7!) 

A gas channel width fraction 

a transfer coefficient for the reaction 

B step size 

y conjugate-gradient coefficient 

X Bruggemann coefficient 

€GDL porosity of gas diffusion layer 

E€V,Ca porosity of catalyst layer 

p electric potential (V) 

(O phase potential (V) 

À membrane water content 
(kmol H20 (kmol SO3~)~!) 

p density of oxygen (kg m~?) 

ogpL electronic conductivity of gas diffusion layer 
(Q7! m7!) 

op electronic conductivity of carbon plate 


(Q7! m7!) 


oscar electronic conductivity of solid catalyst particle 
(Q7! m7!) 

T fluid stress (N m7?) 

‘4 source term in electronic conduction equation 

Subscripts 

a anode 

c cathode/gas channel 

Cat catalyst layer 

eff effective 

GDL gas diffusion layer 

i inlet 

m membrane 

N Nafion 

o outlet 

P carbon plate 

ref reference 

S solid catalyst particles 


From a survey of existing numerical studies, it is readily 
noted that most of the reports were focused on simulation of 
the transport phenomenon in a PEM fuel cell, or at most on 
parametric study of the effects of physical variables. Although 
there is a growing body of literature on the PEM fuel cells, very 
few reports dealt with the optimization of the design parame- 
ters. Among these studies, Mohamed and Jenkins [9] used the 
genetic algorithm to determine the configuration of the PEM 
fuel. Grujicic and Chittajallu [10] proposed a non-linear con- 
strained optimization procedure used in the cathode design in 
order to maximize the average current density at a fixed voltage 
in a PEM fuel cell with interdigitated gas distributors. Optimiza- 
tion was performed for the better catalyst layer composition 
by Secanell et al. [11] using a two-dimensional model and a 
multivariable optimization technique. 

Recently, Lin et al. [12] investigated the optimization of the 
proton exchange membrane fuel cell based on a two-dimensional 
cathode model. Optimization method was employed to deter- 
mine the optimal combination of the design parameters, 
including the channel width ratio, the porosity of GDL and the 
porosity of the catalyst layer. However, the two-dimensional 
half-cell model used by the authors is rather straightforward 
and not able to simulate the complex three-dimensional physical 
transport phenomena in a PEM fuel cell. 

In theory, the task of optimization is basically performed to 
search for a combination of the most appropriate values of a 
number of design parameters so that a defined objective function 
could be minimized. Nowadays, numerous useful optimization 
methods are available, for example, the genetic algorithm [9], the 
conjugate-gradient method (CGM) [13], the Newton-Raphson 
method [14], the adaptive weighting input estimation method 
[15] and so on. Modifications of the traditional conjugate- 
gradient method were proposed by Cheng and Chang [13], and 
the approach is named the simplified conjugate-gradient method 
(SCGM). In the SCGM method, there is no need to solve the 
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simultaneous equations for the step sizes and the complexity 
in mathematical manipulation and the constraint of the form of 
the objective function can then be avoided. Lately, Cheng and 
Chang [16,17] used the concept of the SCGM method to develop 
an inverse method for predicting the internal temperature distri- 
bution of PEM fuel cells. Results show that the outer surface 
temperature data can be used to predict the internal temperature 
distribution at the interface between the carbon plate and the 
MEA. 

Under these circumstances, in the present study, a multi- 
parameter optimization for the geometric parameters of the 
PEM fuel cell based on a three-dimensional full cell model 
is attempted. Schematic of a single-cell PEMFC with parallel 
straight channels on both the anode and the cathode carbon 
plates is shown in Fig. 1. The single-cell PEMFC consists of 
a carbon plate, a gas diffusion layer (GDL), a catalyst layer, for 
each of the anode and the cathode sides, as well as a PEM mem- 
brane at the center. Among these possible influential geometric 
parameters, the gas channel width fraction (A), the gas channel 
height (h) and the thickness of the gas diffusion layer (tgpL) are 
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Fig. 1. Physical model of a single-cell proton exchange membrane fuel cell. (a) 
Single-cell model with parallel straight channels and (b) gas channel geometry. 


(b) 


recognized as the influential factors affecting the behavior of 
the fuel cells. Therefore, these three parameters are selected as 
the design parameters to be optimized in the present study. For 
other physical and geometrical parameters of the cell shown in 
this figure, their values are kept constant throughout the study, 
which are given in Table 1. 

The SCGM method [13] is employed herein as the optimiza- 
tion search scheme. This study is aimed at optimization of these 
three parameters so that the best performance of the PEMFC 
in prescribed operating conditions can be obtained. However, 
it is noted that the present optimization approach is not limited 
to the present group of designed parameters. When necessary, 
more designed parameters may be readily applied. 


2. Optimization approach 


The approach is developed by integrating a direct problem 
solver and an optimizer. The theories of the two solvers are 
briefly described in the following. 


2.1. Direct problem solver 


The following assumptions are made in prior to the construc- 
tion of the governing equations of the three-dimensional full-cell 
model: 


(1) The operating temperature is uniform and constant in the 
solution domain, and in this study, the operating tempera- 
ture and pressure are assigned to be 323 K and 101,325 Pa, 
respectively. 

(2) The fuel cell operates at steady state, and the operating 
voltage is fixed at 0.7 V. 

(3) The gas flows are laminar and compressible in the gas chan- 
nels and the porous layers, and effects of gravity on the flow 
motion are negligible. 

(4) Thermodynamic and electrochemical properties of the gases 
and the solid materials of the fuel cell components are 
assumed constant. 

(5) The gas diffusion layers and the catalyst layers are homo- 
geneous, isotropic porous media. 

(6) Phase change process is ignored and water exists in the entire 
fuel cell only in vapor phase in the simulation. 


2.1.1. Electronic conduction in carbon plates, GDLs and 
catalyst layers 

The carbon plate, the GDL and the catalyst layers all serve 
as conductors for electric current. In the carbon plate and the 
GDL, there is no electrochemical reaction taking place. The 
electrochemical reactions only take place in the catalyst layers. 
Due to the electrochemical reactions activated in the catalyst 
layers, there exists a source term in the electronic conduction 
equation. Therefore, in these three components the electronic 
conduction equation may be expressed in general form as 


VoVo) =t (1) 


where ¢ is the electronic potential and o®f denotes the effec- 
tive electronic conductivity. For the carbon plates, the effective 
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Table 1 

Fixed parameters of base case 

Parameter Symbol Value 
Electronic conductivity of catalyst layer (Q27! m7!) OS.Cat 300 
Effective ionic conductivity of catalyst layer (27! m7!) I; a 4.2 

Effective mass diffusivity of oxygen in catalyst layer (cm? s7!) DE cat 1.9546 x 1073 
Effective mass diffusivity of hydrogen in catalyst layer (cm? s7!) Dp Cat 0.985 x 10-3 
Effective mass diffusivity of oxygen in GDL (cm? s7!) DĚ GDL 1.845 x 107? 
Effective mass diffusivity of hydrogen in GDL (cm? s7!) Dit cpt 0.930 x 107? 
Effective electronic conductivity of catalyst layer (Q7! m7!) o 135.265 
Bruggemann coefficient of gas diffusion layer XGDL 15 
Bruggemann coefficient of catalyst layer XCat 1.5 

GDL permeability (m?) kepi 1.76 x 107!! 
Catalyst layer permeability (m?) Keat 1.76 x 1071! 
Anode inlet gas velocity (m s7!) Va 0.3 

Cathode inlet gas velocity (m s7!) Ve 0.5 

Anode inlet mass fraction of H2 Yup i 0.445 
Anode inlet mass fraction of H20 Yu 0,ai 0.555 
Cathode inlet mass fraction of O2 Yo). 0.212 
Cathode inlet mass fraction of N2 Yn) Ji 0.709 
Cathode inlet mass fraction of H2O Yu 0,ci 0.079 
Operating voltage (V) v 0.7 

Anodic charge transfer coefficients for anode reaction Qaa 0.5 

Cathodic charge transfer coefficients for anode reaction @a,c 0.5 

Anodic charge transfer coefficients for cathode reaction Qca 1.5 

Cathodic charge transfer coefficients for cathode reaction Qcc 1.5 
Reference transfer current density at anode (A m~?) Jaret 9.0 x 108 
Reference transfer current density at cathode (A m~?) Jeet 250 

Anode concentration parameter Ta 1.0 

Cathode concentration parameter Te 1.0 

Porosity of catalyst layer EV,Cat 0.112 
Porosity of gas diffusion layer EGDL 0.5 
Thickness of membrane (m) tm 1.78 x 1074 
Electronic conductivity of gas diffusion layer (Q7! m7!) OGDL 300 
Electronic conductivity of carbon plate (Q7! m7!) op 4000 

Ionic conductivity of Nafion (Q7! m7!) TYN,Cat 25.56 
Operating temperature (K) T 323 

Inlet and outlet pressure (Pa) P;, Po 101,325 

Gas channel length (m) I, 5.0 x 107? 
Width of a module (m) lp 1.0 x 1073 
Thickness of catalyst layer (m) fCat 1.0 x 1075 
Thickness of carbon plate (m) tp 2.0 x 1073 
Height of gas channel (m) h 1.0 x 1073 


electronic conductivity of is treated as a constant (op). For the 
porous gas diffusion and catalyst layers, the magnitude of the 
effective electronic conductivity is determined in terms of the 
porosity in accordance with the Bruggeman’s equation as 


ff 

OGpL = (1 — £anL)*™™oGpL (2a) 
ff 3 

Seat = (1 — Ey, Cat)” 0s, cat (2b) 


It is noticed that the effective electronic conductivity of the 
catalyst layer (of) increases with the volume fraction of solid 
catalyst particles in catalyst layer (€s cat). In addition, the source 
term ¢ in Eq. (1) is zero for the carbon plates or the GDL without 
electrochemical reaction. When the catalyst layers is dealt with, 
¢ is equal to the transfer current densities —ja and je for the anode 
and the cathode, respectively. The transfer current densities at the 
anode and the cathode are calculated based on the Butler- Volmer 


condition as 


; f Yp, f Qaa F —Aac F 
= x | eX ex 
Ja = Ja,ref Fet p RT n p RT n 


(3a) 


. p Yo, ne Ac al Acc F 
= X | eX ex 
Je Je,ref Yozef p RT n p RT n 


(3b) 


where jaref and jc ref are the reference transfer current densities, 
and r, and re the concentration parameters, at the anode and the 
cathode, respectively. In addition, œa, and &a, are the anodic 
and the cathodic charge transfer coefficients for anode reaction, 
and Qe and ac. the anodic and the cathodic charge transfer 
coefficients for cathode reaction. Yo, ref and YH, ref are the ref- 
erence concentrations at which the exchange current densities 
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were obtained. All the reference values of these coefficients are 
already given in Table 1. Note that the reference current densi- 
ties, transfer coefficients and reference concentrations have been 
set so as to fit the experimental data. 


2.1.2. Ionic conduction in catalyst layers and proton 
exchange membrane 

The protons of hydrogen (H*) are generated in the anode 
catalyst layer and then migrate through the proton exchange 
membrane toward the cathode catalyst layer to react with oxygen 
and the returning electrons. Therefore, ionic conduction occurs 
in the catalyst layers and the proton exchange membrane. 

The ionic conduction equation for determination of the phase 
potential (o) in the catalyst layers is 


~VTEV Gcat) = 6 (4) 


where @cat is the phase potential and I~ os is the effective ionic 
conductivity of the catalyst layers. Note that the source term ¢ is 
identical to that of Eq. (1). The effective ionic conductivity I~ bes 
is evaluated by using 


eh, = (1 = evicat)*" I5 Cat (5) 


On the other hand, the ionic conduction equation for the 
membrane is expressed as 


V(I nVm) = 0 (6) 


where 9m is the phase potential and Im is the ionic conductivity 
of the membrane. The expression for ionic conductivity of the 
membrane has been proposed by Springer et al. [18] as: 


1 1 
In(T) = 12 —- = .005239A — 0.0032 
m(T) exp| s (= 7) | (0.005 392. — 0.00326) 
(7) 


where the water content of the membrane (A) is dependent on 
the dimensionless water activity in the membrane (a) according 
to experimental data. The activity is assumed uniform over the 
membrane, which is defined as the gas partial pressure of water 
at equilibrium with the membrane divided by the saturation 
pressure of water vapor at the operating temperature. That is 


Xm oP 
a = z 
Prat 


(8) 


where the saturation pressure of water vapor can be computed 
from 
logio Psat = —2.1794 + 0.02953(T — 273.15) — 9.1837 
x 10->(T — 273.15)? + 1.4454 
M10 T= 273.15) (9) 
with T in K. The correlation between A and a is derived as 
© J 0.043 + 17.81a — 39.85a* + 36a? O<a<1 a0) 
~ | 14+1.4a-— 1) l<a<3 


Based on the solutions for the phase potential (g), the distribu- 
tion of local ionic current density (i;) may be further determined 


with 


i=-—IVg (11) 


2.1.3. Mass, momentum and species transport in gas 
channels, GDLs and catalyst layers 

The governing equations for the gas flows in the gas channels 
and the porous layers are the conservation equations of mass, 
momentum and species. Application of the flow channels and 
the gas diffusion layers in the fuel cells is to allow the reactant 
gases to diffuse uniformly into the catalyst layers but prevent 
excessive increase in electrical resistance against the released 
electrons from the catalyst layers. The porous layers include the 
gas diffusion layers and the catalyst layers, in which the mass, 
momentum and species equations for the flows are derived based 
on the non-Darcy law. These equations are given below: 


Mass equation: 


V(epU) = 0 (12) 
Momentum equation: 
> 2 
— U 
V(eoUU) = —eV P + V(et) + £) 7 (13) 
Species equation: 
V(epUY;) = V(D$ËVY;) + Si (14) 


where £ and k are the porosity and the permeability of the porous 
layers, respectively; t represents the fluid stress tensor. Fuel gas 
at the anode is considered to contain Hz and H20, and oxidant 
gas at the cathode contains O2, N2 and H20. Y; denotes the 
species mass fraction of species i. Note that for the pure fluid 
flows in the gas channels, the values of ¢ and k are given with 
£=] and k=00. 

Calculation for the effective mass diffusivity of species i 
(Dstt) is based on Bruggeman’s equation, that is 


DS" = (e)¥ D; (15) 


where D; is the mass diffusivity of species i and x is the Brugge- 
mann coefficient of the porous layer. 


2.1.4. Electrochemical reaction in catalyst layers 

In Eq. (14), the term S; stands for the source term of species i 
coming from the electrochemical reactions in the catalyst layers. 
Note that H2 and O2 are consumed by the reactions in the catalyst 
layers of anode and cathode, respectively, and H20 is the product 
of the reaction at the cathode. Therefore, the source terms for 
different species are 


Si = 55 (16a) 
j 
So. = Tp (16b) 
j 
Smo = E (16c) 


where ja and ję are calculated with Eqs. (3a) and (3b). 
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Owing to the flexibility that acommercial computational fluid 
dynamics code may have in dealing with the complex three- 
dimensional problems, these above set of conservation equations 
along with the property equations are solved by adopting a finite- 
volume scheme on structured grids within the framework of the 
commercial code CFD-ACE + released by ESI-CFD Inc. 


2.1.5. Boundary conditions 

For the pattern of straight parallel channels shown in Fig. 1, 
a periodic transport phenomenon is expected to develop from 
channel to channel as long as the number of gas channels is large 
enough. Thus, one only needs to deal with the solution domain 
of a single module, with the periodic boundary conditions on the 
faces indicated by the dashed lines. In this study, the width of a 
module (/p) is set to be 1 mm. Therefore, for the faces of solution 
domain at x =0 and 1 mm, a periodic boundary condition for Ü, 
Y, ¢ and ¢ is prescribed. 

As to the inlet conditions, the inlet velocities for the anode gas 
(hydrogen) and the cathode gas (air) in the gas channel, Va and 
Ve, are fixed at 0.3 and 0.5ms~!, respectively. The inlet mass 
fractions for all the gas species are specified, and their values 
are provided in Table 1. On the other hand, for the boundary 
conditions at the exit, it is assumed that the flow becomes fully 
developed downstream. Thus, the exit boundary conditions may 
be described as 


Y; 
ape Y W y a7) 


Solid walls of the gas channels or the carbon plates are imper- 
meable, and hence, on the solid walls, dY;/dn = 0 and Ù = 0. In 
addition, since both the gas channels and the membrane are insu- 
lators against electrons; at the boundaries of the gas channels 
and the membrane, the normal gradient of the electric poten- 
tial is zero, that is, d¢/dn = 0. Furthermore, the ionic conduction 
takes place only inside the catalyst layer and proton exchange 
membrane, and the protons are not allowed to transfer into other 
components; therefore, on the faces of the catalyst layer and the 
membrane the normal gradient of phase potential is assigned to 
be zero. 

The average current density (I) of fuel cells is calculated by 
integration of local electric current density distribution on the 
carbon plate surface at z = 0 divided by the carbon plate surface 
area as 


1 A 
i=; f i dA (18) 


where i, is the local electric current density on the carbon plate 
surface and A is the carbon plate surface area. In the mean time, 
the gas channel width fraction (A) is defined by the ratio of the 
gas channel width (/c) to the width of a module (/p), that is 


=| (19) 


It is important to note that the direct solver used in the present 
work is a single-phase fuel cell model with electron transport 
explicitly taken into account. The inclusion of electron trans- 
port in the fuel cell modeling is required particularly when the 


effect of the gas channel width fraction is considered. For further 
information about the modeling for the electron transport, one 
may refer to the work of Meng and Wang [19]. 


2.2. Optimizer 


In this study, the SCGM method [10] is employed to con- 
struct the optimizer for seeking the optimal values of the design 
parameters. The objective function in conjunction with the opti- 
mization process is defined in the following: 


J(ex) =1/P, k=1-3 (20) 


where eg (k=1-3) represent the design parameters, i.e. gas 
channel width fraction (A), the gas channel height (h) and the 
thickness of the gas diffusion layer (t¢pL). The power density of 
the fuel cell is the product of the cell voltage (V) and the electric 
current density (I) as P=I x V. In this manner, the power den- 
sity (P) is expected to reach a maximum as the magnitude of the 
objective function (J) is minimized. In the SCGM method, the 
values of the step size (xz) during the searching procedure are 
fixed. The initial guess for each design parameter is made first, 
and in the successive steps the conjugate-gradient coefficients 
and the searching directions are evaluated. The new design vari- 
ables are continuously updated. When the objective function 
reaches a minimum, the optimization process is completed and 
the searching procedure is then terminated. 
The procedure in the optimizer is presented as follows: 


(1) Make the initial guess for the designed parameters eg, assign 
the values to the step sizes Bx. 

(2) Specify all boundary conditions, and then solve the physical 
problem by the numerical simulator. 

(3) Calculate the objective function J (eg, k = 1-3). If the conver- 
gence criterion is satisfied as the objective function reaches 
a minimum in the iterating process, the solution process is 
terminated. Otherwise, proceed to step (4). 


(4) Calculate the gradient functions of the function J by the 
sensitivity analysis of numerical direct differentiation [13]: 
ay” AJ” 

= , k=123 (21) 
dek Aek 

(5) Calculate the conjugate-gradient coefficients y¢ (k = 1-3): 

dJ/dex)" 7? 
m | k=1-3 (22) 
(8J/ 3er)” 
For the first step with n=0, y? =0. 
(6) Calculate the searching directions &/(k = 1-3): 
ðJ” 
ET! = + Ke, k= 13 (23) 
dek 

(7) Update new design variables by 

etl e — pert! k=1-3 (24) 


and then return to step (2). 
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Optimizer 
(SCGM method) 
(self-developed) 


Direct Problem Solver 
( CFD-ACE+) 


Pre-Processor 
(CFD-GEOM) 


Fig. 2. Connection between optimizer and direct problem solver. 


2.3. Connection between direct problem solver and 
optimizer 


The self-developed optimizer and the computational fluid 
dynamics code with the pre-processor are connected through 
an interface program. The interface programs are developed by 
using the Python language. Connection among the optimizer 
and the direct problem solver is shown in Fig. 2. Through the 
Python-interface connection, the message of any change in the 
design parameters suggested by the optimizer is firstly sent to 
the pre-processor for updating the grid system and the bound- 
ary conditions. Secondly, the CFD code is executed based on 
the updated grid system and boundary conditions to yield the 
numerical solutions for distribution of U , Y, @ and ¢ in all the 
components of the fuel cell. The solutions obtained from the 
direct problem solver are then used to calculate the value of the 
objective function, which is further transmitted back to the opti- 
mizer for calculating the consecutive searching directions. Flow 
chart of the optimization process is plotted in Fig. 3. 


2.4. Experiments for validation 


To verify the numerical predictions of the direct problem 
solver, a single-cell fuel cell module with cross-sectional area 
of 22cm x 22cm and active area of 14cm x 14cm has been 
installed and tested. Two 3-mm thick composite graphite plates 
are used as the flow fields and current conductors. A five- 
layer membrane exchange assembly (MEA) is placed at the 
center between the two composite graphite plates. The mem- 
brane used in this MEA is Nafion® 112 and the gas diffusion 
layers are made of carbon fiber paper. Meanwhile, in the exper- 
iments the volumetric flow rates of the hydrogen and the air are 
typically maintained at 1500 and 2500 sccm, corresponding to 
V,=0.3ms—! and V.=0.5ms~!, respectively. The hydrogen 
and the air are humidified to reach 100% relative humidity and 
80°C temperature before entering the gas flow channels. The 
operating and geometrical conditions given to the direct prob- 
lem solver are consistent with the experiments. The fuel cell 
temperature is maintained at 343 K. 

Fig. 4 shows a comparison between the numerical predictions 
and the experiment data for the polarization curves of the single- 
cell fuel cell module. It is observed that the numerical predictions 
closely agree with the experiments, especially in the regime of 
1<3800A m~°. Over this regime, the experimental data appear 
to experience a greater drop in the operating voltage than the pre- 
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Fig. 3. Flow chart of the optimization process. 
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Fig. 4. Comparison in polarization curves between predictions and experiments. 
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dicted. The greater voltage drop is attributed to the concentration 
loss due to insufficient supply of O2 and H2 in experiments at the 
cathode and the anode, respectively. Nevertheless, in this figure 
it is found that at V=0.7 V the numerical and the experimental 
data for the current density are in close agreement (at approx- 
imately 2500 A m°). Thus, the typical value of the operating 
voltage of the fuel cell is fixed at 0.7 V in the present study to 
determine the optimal parameters. 


3. Results and discussion 
3.1. Parametric study 
In addition to the experimental demonstration shown earlier, 


the direct problem solver’s capability of providing numerical 
predictions under various conditions is further studied herein. 


Cathode 


x x 


x=0 mm x=Ilmm 
(a) (b) 


Numerical predictions of the mass fraction distributions of 
oxygen and hydrogen at the cathode and the anode sides, respec- 
tively, on cross section at y= 25 mm with different channel width 
fractions (A =0.6, 0.4 and 0.2) are shown in Fig. 5. In the mean 
time, planar distributions of oxygen mass fraction distribution 
at the center plane of the cathode catalyst layer with different 
channel width fractions are displayed in Fig. 6. The interface 
between the gas channel and the rib is indicated by the dashed 
lines. It is expected that a larger channel width fraction leads to 
sufficient gases supply into the catalyst layer for electrochemi- 
cal reaction. The results shown in the two figures clearly match 
this expectation. As observed in Figs. 5 and 6, the concentra- 
tions of the reactants are higher in the gas channels at A =0.6, 
whereas at A =0.2 the concentrations of the reactants are sig- 
nificantly increased. However, note that a larger channel width 
also gives narrower ribs which serve as electron-conduction pas- 


(c) 


Fig. 5. Mass fraction distributions on cross section at y = 25 mm with different channel width fractions (dashed line indicates that interface between channel and rib): 


(a) A=0.6, (b) A=0.4 and (c) A=0.2. 
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Fig. 6. Planar distributions of oxygen mass fraction distribution at the center 
plane of the cathode catalyst layer with different channel width fractions (dashed 
line indicates that interface between channel and rib): (a) A=0.6, (b) A=0.4 
and (c) A =0.2. 


sages; therefore, an increase in the gas channel fraction is always 
accompanied by an increase in the ohmic impedance. With these 
two opposing effects, there should exist a desired optimal value 
of the gas channel width fraction. 

Fig. 7 conveys the distributions of oxygen mass fraction 
distribution on the center plane of the cathode catalyst layer 
with different thicknesses of gas diffusion layer (tgpL). The 
thickness of the gas diffusion layer is assigned to be 1 x 1074, 
2 x 10-4 and5 x 1074 m. It is observed in Fig. 7(c) that a thicker 
gas diffusion layer (say, 5 x 1074m) may allow the reactant 
gas to diffuse more uniformly in the entire gas diffusion layer 
before penetrating into the catalyst layer. Hence, on the center 
plane of the cathode catalyst layer, a very uniform concentra- 
tion distribution can be seen. On the other hand, in Fig. 7(a) 
for the case at tgp, =1 x 1074m, the oxygen concentration 
is appreciably higher in the area adjacent to the gas channel 
than in the area adjacent to the rib. However, an increase in 
the thicknesses of gas diffusion layer also tends to increase the 
electric resistance against the conduction of electrons. Here, 


| y= 0mm 


x=Omm x=1mm 


(a) (b) (c) 


Fig. 7. Planar distributions of oxygen mass fraction distribution on the center 
plane of the cathode catalyst layer with different thicknesses of gas diffu- 
sion layer (dashed line indicates that interface between channel and rib): (a) 
top. = 1 x 1074 m, (b) tapk =2 x 1074 m and (c) tgp, =5 x 1074 m. 


one finds another influential parameter that is worth further 
study. 

Plotted in Fig. 8 are the local current density distributions 
on the center plane of the cathode catalyst layer with different 
channel heights (A). In this figure, the gas channel height is 
assigned to be 4 x 1074, 8 x 1074 or 16 x 1074m. Mass flow 
rates of the reactant gases can be elevated by increasing the 
height of gas channels, so as to provide more gases for reaction; 
however, higher channels may also slightly elevate the ohmic 
impedance of the fuel cell. It is indeed observed in Fig. 8 that 
the magnitude of the local current density is higher in the area 
adjacent to the rib since the rib serves as passage for electron 
conduction. Meanwhile, local current density becomes stronger 
on the entire plane when the channel height becomes lower. 


3.2. Optimization of PEM fuel cell 


In the parametric studies considered above, only one param- 
eter is changed while others are fixed. In fact, multi-parameter 
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x= | mm 


x= 0 mm 
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Fig. 8. Local current density distributions on the center plane of the cathode 
catalyst layer with different channel heights (dashed line indicates that inter- 
face between channel and rib): (a) h=4 x 1074 m, (b) h=8 x 1074m and (c) 
h=16 x 1074 m. 


optimization is more meaningful in practices. Design of the geo- 
metric parameters is now carried out to optimize the performance 
of PEMFC by integrating the direct problem solver with the 
optimizer. 

Attention is now drawn to the results of optimization for the 
geometric parameters, A, tgp and h. Trajectory of the optimiza- 
tion process seeking an optimal set of these design parameters 
is shown in Fig. 9. The optimal combination is sought in the 
parameter space (A, tGpL, 4). The optimization process is started 
with the original design of A =0.8, tgpL = 1 mm and h = 0.8 mm. 
It can be observed in Fig. 9 that the SCGM optimizer effec- 
tively brings the searching process to approach the optimal state. 
In this case, it is found that the iteration is converged after 
about 30 iterations. The computation time required is approx- 
imately 60h on a personal computer with an AMD-3GMHz 
CPU. Based on the results shown in Fig. 9, it is found that at 
A=0.3925, tgpL = 0.176 mm and A} = 1.2034 mm, the optimiza- 
tion process reaches a maximum power of 1857 W m~?. The 
variation of power density during the optimization process is 
shown in Fig. 10. Note that the power density associated with 
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Fig. 9. Trajectory of the optimization process for a PEM fuel cell seeking an 
optimal set of designed parameters, A, tgp and h. 


the original design is approximately 1400 W m7. The power 
density then approaches its maximum value when the optimal 
parameters are obtained. As a result, the optimization process 
leads to an increase by 457 W m~? (32%) in the power density. 

Fig. 11 shows the comparison in performance curve between 
the original and the optimal design obtained by the optimizer. 
Results show that the present approach can yield an optimal 
design having higher performance than the original one. Based 
on the features of the curves, it is recognized that the optimal 
design is greatly improved not only in the ohmic loss but also in 
the limiting current density. 

In order to ensure the robustness of the optimizer, Fig. 12 
shows the convergence of the optimization process from three 
different initial combinations of the three parameters, A, tcp. 
and h. It can be observed in this figure that the present approach 
leads to the same optimal set of parameters from different initial 
guesses. 

Note that the suggested optimal parameters are obtained 
under the base case conditions considered in the present study. 
When the conditions are changed, the optimal values of these 
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Fig. 10. Variation in power density during the optimization process. 
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Fig. 11. Comparison in the polarization curves between the original and the 
optimal designs. 


Initial guess 1 


0.4 ‘cp. [mm] 


Fig. 12. Trajectories of the optimization processes from different sets of initial 
guesses for A, tgpL and h. 


parameters may be altered. It is also important to note that exten- 
sion of the present method to the design with more parameters 
is a worthwhile work in the future, and a more comprehensive 
model would also be desired. 


4. Conclusions 


The present study is concerned with the development of an 
numerical approach which is applied for optimization of the 
geometric parameters of the PEM fuel cells. The approach is suc- 
cessfully developed by integrating a direct problem solver with 
an optimizer. A commercial computational fluid dynamics code 
is used as the direct problem solver, which is used to simulate 
the three-dimensional mass, momentum and species transport 


phenomena as well as the electron- and proton-transfer process 
taking place in a PEMFC. On the other hand, the SCGM method 
[10] is employed to build the optimizer, which is combined with 
the direct problem solver in order to seek the optimal geometric 
parameters. In this study, geometric parameters including the 
gas channel width fraction, the gas channel height and the thick- 
ness of the gas diffusion layer are regarded as the major factors 
to optimize. 

A 22cm x 22cm single-cell fuel cell module with 14cm x 
14cm active area has been installed and tested. The numerical 
predictions by the direct problem solver have been verified by 
the experiments. 

For the particular conditions considered in this study, 
it is found that that at A=0.3925, tgp, =0.176mm and 
h=1.2034mm, the optimization process reaches a maximum 
power density of 1857 Wm’, with an increase of 457 W m~? 
(32%) compared with the original design. 

It is also found that the present approach can be applied 
to determine the optimal set of geometric parameters, and the 
search process is robust and always leads to a unique final solu- 
tion regardless of the initial guess. 

Note that the suggested optimal parameters are obtained 
under the base case conditions considered in the present study. 
When the conditions are changed, the optimal values of these 
parameters may be altered. 


Acknowledgements 


The authors would like to thank the National Science Coun- 
cil, Taiwan, ROC, for their financial support under Grant: NSC 
93-2212-E-036-001. Authors would also like to thank Tatung 
Company for their support of the experimental system. 


References 


[1] J.T. Wang, R.F. Savinell, Electrochim. Acta 37 (1992) 2737-2745. 
[2] J.S. Yi, T.V. Nguyen, J. Electrochem. Soc. 146 (1999) 38-45. 
[3] V. Gurau, F. Barbir, H. Liu, J. Electrochem. Soc. 147 (2000) 4485-4493. 
[4] S. Um, C.Y. Wang, K.S. Chen, J. Electrochem. Soc. 147 (2000) 4485-4493. 
[5] W. Sun, B.A. Peppley, K. Karana, J. Power Sources 144 (2005) 42-53. 
[6] E. Hontanon, M.J. Escudero, C. Bautista, P.L. Garcia-Ybarra, L. Daza, J. 
Power Sources 86 (2000) 363-368. 
[7] T. Berning, N. Djilali, J. Power Sources 124 (2003) 440-452. 
[8] C.Y. Wang, Chem. Rev. 104 (2004) 4727-4766. 
[9] I. Mohamed, N. Jenkins, J. Power Sources 131 (2004) 142-146. 
[10] M. Grujicic, K.M. Chittajallu, Chem. Eng. Sci. 59 (2004) 5883-5895. 
[11] M. Secanell, et al., Electrochim. Acta 52 (2007) 2668-2682. 
[12] H.H. Lin, C.H. Cheng, C.Y. Soong, F. Chen, W.M. Yan, J. Power Sources 
162 (2006) 246-254. 
[13] C.H. Cheng, M.H. Chang, Numer. Heat Transfer B 43 (2003) 489-507. 
[14] Y. Hung, N.G. Zamani, Appl. Math. Comput. 80 (1996) 155-179. 
[15] P.C. Tuan, M.C. Ju, Numer. Heat Transfer B 37 (2000) 247-265. 
[16] C.H. Cheng, M.H. Chang, J. Power Sources 139 (2005) 115-125. 
[17] M.H. Chang, C.H. Cheng, J. Power Sources 142 (2005) 200-210. 
[18] T.E. Springer, T.A. Zawodzinski, S. Gottesfeld, J. Electrochem. Soc. 138 
(1991) 2334-2342. 
[19] H. Meng, C.Y. Wang, J. Electrochem. Soc. 151 (2004) A358-A367. 


